! metropolis MC of 2D xy-model sixhoney  lattics 
! metropolis ising+XY model on the square lattice   2021/9/3
! 把XY 模型 超流密度 怎么测量，怎么定义，什么样
! 平衡以后，相变点之前,isp()  xisp() 
!    + - + -   isp (i)   
!    - + + —
!    - + - +
!    - - - -
!    pi 0.5*PI  0.6*PI  1*PI 
!
!
!
!
      implicit real*8 (a-h,o-z)                                         
      real*8 M_2,M_4
      save

       read*,n1,ntoss,iseed,ncyc,nm,temp,Jc,Q
c      print*,n1,ntoss,iseed,ncyc,nm,temp,Jc,Q 
c      Tc=0.8935
c      print*,n1,ntoss,iseed,ncyc,nm,temp
c      if ((n1.eq.0).or.(n1.gt.80).or.(iseed.eq.0)) stop
c      initialize, equilibrate
12    continue
      call initial(n1,iseed)! lattice, nbor, xisp, isp 
!       call mcmxywolff(delta,Q,ntoss*nm,temp)
      call mcmising_xymetroflip(Jc,q,nm*ntoss,temp)
      arm=0              ! average relative magnetization
      asm=0              ! average square magnetization
      ave=0              ! average relative energy
      ase=0              ! average square energy
      M_2=0              !update9.15
      M_4=0              !update9.15
c     simulate and sample
      do 11 ip=1,ncyc                                                   
!       call mcmxywolff(delta,Q,nm,temp)
       call mcmising_xymetroflip(Jc,q,nm,temp)
       call sampleGxy(Jc,Q,rem,rem_2,enr)

c       call pry(temp,Tc)
c       call prx(n1) 
      sqm=rem*rem            ! calculate square magnetization
      arm=arm+rem            ! accumulate relative magnetization
      asm=asm+sqm            ! accumulate square magnetization
      ave=ave+enr            ! accumulate energy density
      ase=ase+enr*enr        ! accumulate square energy
      M_2=M_2+rem_2*rem_2    ! update 9.15
      M_4=M_4+rem_2**4       ! update 9.15
 11   continue                                                          
c normalize, analyze, and print output
      arm=arm/ncyc
      asm=asm/ncyc
      ave=ave/ncyc
      ase=ase/ncyc
      M_2=M_2/ncyc           !update9.15
      M_4=M_4/ncyc           !update9.15
      susp=asm-arm*arm
      speh=(ase-ave*ave)/temp/temp  
      tbr=(M_2**2)/M_4       !update9.15    
      open(2,file='res.dat',status='unknown',access='append') 
      print'(8f12.5)',temp,arm,asm,susp,ase,ave,speh,tbr
	  write(2, 13) temp,arm,asm,susp,ase,ave,speh,tbr
      close(2) 
 13   format(8f10.6)
c      read(1,*),temp
      read *, temp 
c     print*,'temp=',temp 
      if (temp.ne.-1.d0) goto 12       !stop
      close(1) 
      end

         
      subroutine prx(m1)
      implicit real*8(a-h,o-z)
      common/xyc/xisp(6400) 
      dimension C(6400)
      C(1:6400)=0     
      nsq=m1*m1
      
      open(1,file='8xtrain.dat',access='append') 

      do 111 ispin = 1,nsq        
c          write(1,*) cos(xisp(ispin)), sin(xisp(ispin))  
          C(2*ispin-1)=cos(xisp(ispin))
          C(2*ispin)= sin(xisp(ispin)) 
111   continue
      write(1,35) C(1:2*nsq)  
35    format(400f10.5)           
      close(1)  
      end


      subroutine pry(temp,Tc)
c subroutine for y.dat
      implicit real*8 (a-h,o-z) 
      open(1,file= '8ytrain.dat',access='append')  
      if (temp<Tc) then 
          write(1,*) 0     
      else
          write(1,*) 1 
      endif
      close(1)
      end
      
      subroutine mcmxymetroflip(delta,q,nsteps,temp)                   
c subroutine for metropolis sweeps
      implicit real*8 (a-h,o-z)
      common/lattice/ isp(6400),nbor(4,6400),n1,nsq,nsqh,msign,nwc,ngc
      common /lookup/ ppr(-6:6)
chttp://www.lps.ens.fr/~krauth/images/7/72/Stage_Mayer_Johannes_2015.pdf      
      dimension nstack(6400), kflag(6400)
      parameter (pi=3.1415926d0)
      common/xyc/xisp(6400)
      do 400 isteps=1,nsteps
      do 400 ispin=1,nsq

      ns=rn()*nsq+1
       delta1 = delta
c      print*,delta1,'xxxxxxx' 
      alpha=rn()*2*pi
c      print*,'alpha=',alpha
      
c      dt=0.1*rn()*(2*int(2*rn())-1)!

      Eold= delta1 * (-cos(xisp(ns)-xisp(nbor(4,ns)))) +(1 - delta1)*
     &(-cos(q*xisp(ns)-q*xisp(nbor(4,ns))))  + delta1 * (-cos(xisp(ns)-
     &xisp(nbor(1,ns))))+(1 - delta1)*(-cos(q*xisp(ns)-
     &q*xisp(nbor(1,ns))))+ delta1 * (-cos(xisp(ns)-xisp(nbor(2,ns))))
     &+(1 - delta1)*(-cos(q*xisp(ns)-q*xisp(nbor(2,ns))))+
     &(-cos(xisp(ns)-xisp(nbor(3,ns)))) * delta1 + (1 - delta1)*
     &(-cos(q*xisp(ns)-q*xisp(nbor(3,ns))))

           temp1 =2*alpha- xisp(ns)

      if (temp1.gt.2*pi) then
          temp1=temp1-2*pi
      endif
      if (temp1.lt.0) then
          temp1=temp1+2*pi
      endif

      Enew= delta1*(-cos(temp1-xisp(nbor(4,ns))))+ (1 - delta1)*
     &(-cos(q*temp1-q*xisp(nbor(4,ns))))+delta1*(-cos(temp1-
     &xisp(nbor(1,ns))))+(1 - delta1)*(-cos(q*temp1-q*xisp(nbor(1,ns))))
     &+delta1*(-cos(temp1-xisp(nbor(2,ns))))+(1 -delta1)*(-cos(q*temp1
     &-q*xisp(nbor(2,ns))))+delta1*(-cos(temp1-xisp(nbor(3,ns))))+
     &(1-delta1)*(-cos(q*temp1-q*xisp(nbor(3,ns))))

       deltaE =Enew-Eold



      prb = exp(-deltaE/temp)
      if (rn().lt.prb) then
           xisp(ns)=temp1
      end if






c      q= 1
400   continue
      return  
      end 

      subroutine mcmising_xymetroflip(Jc,q,nsteps,temp)                   
      implicit real*8 (a-h,o-z)
      common/lattice/ isp(6400),nbor(4,6400),n1,nsq,nsqh,msign,nwc,ngc
      common /lookup/ ppr(-6:6)
c     update 2021.9.8    
      dimension nstack(6400), kflag(6400)
      parameter (pi=3.1415926d0)
      common/xyc/xisp(6400)
      do 400 isteps=1,nsteps
      do 400 ispin=1,nsq

      ns=rn()*nsq+1
	  !16*0.1+1=2  1 2 3 4 
	  !            5 6 7 8 
	  !16*0.1+1=2  9 101112 
	  !            13141516
      
      alpha=rn()*2*pi
	  !     2
	  !  5  6  7 
	  !     10
	  ! eold_ij= -Jcouping*(1+ sigma_i * sigma_j) *( cos (theta_i-theta_j))
	  
      t1= -Jc*cos(xisp(ns) - xisp(nbor(4,ns)))
	  t2= (1 + isp(ns)*isp(nbor(4,ns))) 
	  t3=    -Jc*cos(xisp(ns) - xisp(nbor(3,ns)))
	  t4=     (1 + isp(ns)*isp(nbor(3,ns)))
	  t5=    -Jc*cos(xisp(ns) - xisp(nbor(2,ns)))
	  t6=     (1 + isp(ns)*isp(nbor(2,ns)))
	  t7=   -Jc*cos(xisp(ns) - xisp(nbor(1,ns)))
	  t8=    (1 + isp(ns)*isp(nbor(1,ns)))
      Eold= t1*t2+ t3*t4+ t5*t6+ t7*t8 

           temp1 =2*alpha- xisp(ns)

      if (temp1.gt.2*pi) then
          temp1=temp1-2*pi
      endif
      if (temp1.lt.0) then
          temp1=temp1+2*pi
      endif
	  t1= -Jc*cos(temp1 - xisp(nbor(4,ns)))
	  t2= (1 + isp(ns)*isp(nbor(4,ns))) 
	  t3=    -Jc*cos(temp1 - xisp(nbor(3,ns)))
	  t4=     (1 + isp(ns)*isp(nbor(3,ns)))
	  t5=    -Jc*cos(temp1 - xisp(nbor(2,ns)))
	  t6=     (1 + isp(ns)*isp(nbor(2,ns)))
	  t7=   -Jc*cos(temp1 - xisp(nbor(1,ns)))
	  t8=    (1 + isp(ns)*isp(nbor(1,ns)))
      Enew= t1*t2+ t3*t4+ t5*t6+ t7*t8 
      deltaE =Enew-Eold
      prb = exp(-deltaE/temp)
      if (rn().lt.prb) then
           xisp(ns)=temp1
      end if


      ns=rn()*nsq+1
	  !16*0.1+1=2  1 2 3 4 
	  !            5 6 7 8 
	  !16*0.1+1=2  9 101112 
	  !            13141516
      
      
	  !     2
	  !  5  6  7 
	  !     10
	  ! eold_ij= -Jcouping*(1+ sigma_i * sigma_j) *( cos (theta_i-theta_j))
	  
      t1= -Jc*cos(xisp(ns) - xisp(nbor(4,ns)))
	  t2= (1 + isp(ns)*isp(nbor(4,ns))) 
	  t3=    -Jc*cos(xisp(ns) - xisp(nbor(3,ns)))
	  t4=     (1 + isp(ns)*isp(nbor(3,ns)))
	  t5=    -Jc*cos(xisp(ns) - xisp(nbor(2,ns)))
	  t6=     (1 + isp(ns)*isp(nbor(2,ns)))
	  t7=   -Jc*cos(xisp(ns) - xisp(nbor(1,ns)))
	  t8=    (1 + isp(ns)*isp(nbor(1,ns)))
      Eold= t1*t2+ t3*t4+ t5*t6+ t7*t8 

          
      
     
	  t1= -Jc*cos(xisp(ns) - xisp(nbor(4,ns)))
	  t2= (1 - isp(ns)*isp(nbor(4,ns))) 
	  t3=    -Jc*cos(xisp(ns) - xisp(nbor(3,ns)))
	  t4= (1 - isp(ns)*isp(nbor(3,ns)))
	  t5=    -Jc*cos(xisp(ns) - xisp(nbor(2,ns)))
	  t6= (1 - isp(ns)*isp(nbor(2,ns)))
	  t7=   -Jc*cos(xisp(ns) - xisp(nbor(1,ns)))
	  t8= (1 - isp(ns)*isp(nbor(1,ns)))
      Enew= t1*t2+ t3*t4+ t5*t6+ t7*t8 
      deltaE =Enew-Eold
      prb = exp(-deltaE/temp)
      if (rn().lt.prb) then
           isp(ns) = -isp(ns) 
      end if






c      q=1
400   continue
      return  
      end 

      
      subroutine mcmxywolff(delta,q,nsteps,temp)                       
c subroutine for metropolis sweeps
      implicit real*8 (a-h,o-z)                                         
      common/lattice/ isp(6400),nbor(4,6400),n1,nsq,nsqh,msign,nwc,ngc
      common /lookup/ ppr(-6:6)
c     http://www.lps.ens.fr/~krauth/images/7/72/Stage_Mayer_Johannes_2015.pdf      
      dimension nstack(6400), kflag(6400)
      parameter (pi=3.1415926d0)
      common/xyc/xisp(6400)
	  dimension xold(6400)

c    1. Choosing randomly one site k of the lattice, we will call it the seed.
c    2. Choosing randomly an angle α∈ [0,π).
c    3. Flipping the spin at the site k of the lattice around the axis α , that means
c    θ∗ k = 2 α − θ_k.
c    4. Flip each unflipped neighbour i of the site k with probability P(k,i)
c     and repeat this procedure for the site i.

     
      do 200 isteps=1,nsteps
          do k=1,nsq
            kflag(k)=-2 ! -1: unscanned, 1 scanned
			xold(k)=xisp(k)
          enddo       
          isd=0;  k=rn()*nsq+1 ; delta1=delta;  alpha=pi*rn()
          kflag(k)=2 ;    isd=isd+1; nstack(isd)=k;
          xisp(k)=2*alpha-xisp(k)        
          do while (isd.ge.1)
             k=nstack(isd);   isd=isd-1;
             do i=1,4
             if (kflag(nbor(i,k)).eq.-2) then	
!             t1= -Jc*cos(xold(k) - xisp(nbor(i,k)))
!	         t2= (1 + isp(k)*isp(nbor(i,k))) 
!             Eold= t1*t2
!			 t1= -Jc*cos(xisp(k) - xisp(nbor(i,k)))
!	         t2= (1 + isp(k)*isp(nbor(i,k))) 
!             Enew= t1*t2
			 
            Eold= delta1 * (-cos(xold(k)-xisp(nbor(i,k))))
     &+(1 - delta1)*(-cos(q*xold(k)-q*xisp(nbor(i,k))))
             Enew= delta1 * (-cos(xisp(k)-xisp(nbor(i,k))))
     &+(1 - delta1)*(-cos(q*xisp(k)-q*xisp(nbor(i,k))))
                twobetasidotrsjdotr= (Eold-Enew)/temp
          padd=1-min(1.0,exp(twobetasidotrsjdotr))
      
		    if (rn()<padd) then
            isd=isd+1;
            nstack(isd)=nbor(i,k);
            kflag(nbor(i,k))=2;
            xisp(nbor(i,k))=2*alpha-xisp(nbor(i,k))
          
                          if (xisp(nbor(i,k)).gt.2*pi) then
                 xisp(nbor(i,k))=xisp(nbor(i,k))-2*pi
                endif
               if (xisp(nbor(i,k)).lt.0) then
                 xisp(nbor(i,k))=xisp(nbor(i,k))+2*pi
                endif
           
            endif 
		 endif
            enddo  
           id=id+1  
          enddo 
         
100   continue     
       
      
   
        
 200  continue
    
      return                                                            
      end                                                               






      


      subroutine  ihoney(m1,iseed)
      implicit real*8(a-h,o-z) 
      save
      common/xyc/xisp(6400) ! newadded 
      parameter(pi=3.1415926d0) !newadded 
      common /lookup/ ppr(-6:6)
      common/lattice/ isp(6400),nbor(4,6400),n1,nsq,nsqh,msign,nwc,ngc
      nx=m1*2
     
      ny=m1
      nsq=nx*ny
c   define neighbors of each spin
      do 410 ispin=1,nsq
      iy=(ispin-1)/nx+1
      ix=ispin-(iy-1)*nx
      ixp=ix+1-(ix/nx)*nx
      iyp=iy+1-(iy/ny)*ny
      ixm=ix-1+((nx-ix+1)/nx)*nx
      iym=iy-1+((ny-iy+1)/ny)*ny
c      print*,ix,iy,ixm,ixp,iym,iyp
      nbor(1,ispin)=(iy-1)*nx+ixm
      nbor(2,ispin)=(iy-1)*nx+ixp
      if (mod(ispin,2).eq.1) then
      nbor(3,ispin)=(iyp-1)*nx+ixp
      else
      nbor(3,ispin)=(iym-1)*nx+ixm
      endif
410   continue

      if(iseed.eq.0) goto 420        

c initialize random generator
      call ransi(iabs(iseed)+1)
      is=1
     
      do 412 ns=1,nsq
      isp(ns)=is
c      xisp(ns)=rn()*pi*(2*int(2*rn())-1)
      xisp(ns)=0 
412   is=is
c initialize geometric cluster data
      maxi=2**30-1+2**30             ! largest 32-bit integer+
      nrcl=maxi
420   continue

      return
      end

c      print*,'final'
c      do is=1,nsq
c       print*,is,nbor(1:3,is)

       
      subroutine  initial(m1,iseed)
c     initialize for Ising mc                                           
      implicit real*8(a-h,o-z)                                          
      save
      common /lookup/ ppr(-6:6)
      common/lattice/ isp(6400),nbor(4,6400),n1,nsq,nsqh,msign,nwc,ngc
      common/xyc/xisp(6400) ! newadded configurations of the xy model 
      parameter(pi=3.1415926d0) !newadded 
      n1=m1                                                             
      nsq=n1*n1                                                        
      nsqh=nsq/2
      msign=1
c define neighbors of each spin
      do 210 ispin=1,nsq                                                
      iy=(ispin-1)/n1+1                                                
      ix=ispin-(iy-1)*n1                                     
      ixp=ix+1-(ix/n1)*n1                                               
      iyp=iy+1-(iy/n1)*n1                                               
      ixm=ix-1+((n1-ix+1)/n1)*n1                                        
      iym=iy-1+((n1-iy+1)/n1)*n1                                        
      nbor(1,ispin)=(iy -1)*n1+ixm 
      nbor(4,ispin)=(iy -1)*n1+ixp 
      nbor(2,ispin)=(iym-1)*n1+ix  
      nbor(3,ispin)=(iyp-1)*n1+ix
210   continue                                                          
c fill lookup table for local updates and define bond probabilities
      bfm=dexp(1.0d0/temp) 
      bf=1.0d0                                                       
      bfi=bf
      bp=1.0d0-1.0d0/(bfm*bfm)
      gp=1.0d0-1.0d0/(bfm*bfm*bfm*bfm)
      ppr(0)=0.5d0 
      do 220 i=1,6                                                   
      bf=bf*bfm                                                         
      bfi=bfi/bfm                                                       
      ppr(i)=bf/(bf+bfi)
      ppr(-i)=bfi/(bf+bfi)
220   continue                                                          
      pb=1.0d0-1.0d0/(bfm*bfm)
      if(iseed.eq.0) goto 320           
                                   
c initialize random generator
      call ransi(iabs(iseed)+1)
      is=1
      do 112 ns=1,nsq
      isp(ns)=is
      xisp(ns)=2*pi*rn() 
112   is=is
c initialize geometric cluster data
      maxi=2**30-1+2**30             ! largest 32-bit integer+
      nrcl=maxi
320   continue                                                          
      
      return                                                            
      end                                                               


      subroutine sampleGXY(Jc,q,rem,rem_2,enr)
c     compute observables from given spin configuration                 
      implicit real*8 (a-h,o-z)
      common/lattice/ isp(6400),nbor(4,6400),n1,nsq,nsqh,msign,nwc,ngc
      common/xyc/xisp(6400)
      dimension xm(10,3) 
      save
	  xbs = 0 
	  m1  = 0 
	  do ns=1,nsq 
          t1= -Jc*cos(xisp(ns) - xisp(nbor(4,ns)))
	  t2= (1 + isp(ns)*isp(nbor(4,ns))) 
	  t3=    -Jc*cos(xisp(ns) - xisp(nbor(3,ns)))
	  t4=     (1 + isp(ns)*isp(nbor(3,ns)))
	  t5=    -Jc*cos(xisp(ns) - xisp(nbor(2,ns)))
	  t6=     (1 + isp(ns)*isp(nbor(2,ns)))
	  t7=   -Jc*cos(xisp(ns) - xisp(nbor(1,ns)))
	  t8=    (1 + isp(ns)*isp(nbor(1,ns)))
      Eold= t1*t2+ t3*t4+ t5*t6+ t7*t8 
	  xbs = xbs + Eold 
	  m1  = m1 + isp(ns) 
      enddo 
	  xbs= xbs/2
	  rem=abs(float(m1)/nsq)! absolute magnetization per site 
          rem_2=float(m1)/nsq   !update9.15求m
	  !1/2=0  1.0/2 =0.5 
       
      enr=xbs/nsq! average energy per site 
	  
       end 

      subroutine sample(rem,enr)
c     compute observables from given spin configuration                 
      implicit real*8 (a-h,o-z)                                         
      common/lattice/ isp(6400),nbor(4,6400),n1,nsq,nsqh,msign,nwc,ngc
      common/xyc/xisp(6400)
      save
c     sample magnetization 
      m1=0
      m2=0
      do 20 ns=1,nsqh
      m1=m1+isp(ns)
      m2=m2+isp(ns+nsqh)
 20   continue                                                          
c     sample nearest-neighbour sum
      nbs=0
      do 30 ns=1,nsq                                                    
      nbs=nbs+isp(ns)*(isp(nbor(1,ns))+isp(nbor(2,ns)))
 30   continue                                                          

      bsn=nbs
      bsn=bsn/(n1*n1)
      rem=m1+msign*m2
      rem=dabs(rem/(n1*n1))
      enr=bsn
      return                                                            
      end                                                               

      subroutine ransi(iseed)                                 
c initialize shift register random generator 
      implicit real*8(a-h,o-z)            
      save
      parameter (mult=32781,lenr=9689,ifdb=471)
      common/ransrb/ irs(lenr),next(lenr),ipoint,ipoinf
      k=3**18+2*iseed     
      do 100 i=1,lenr      
      k=k*mult               
      irs(i)=k+i/3+i/17       
 100  continue                 
      do 101 i=1,lenr           
 101  next(i)=i+1                
      next(lenr)=1                
      ipoint=1                     
      ipoinf=ifdb+1                
      return                         
      end

      function rn()
c     calculate random number      
      implicit real*8(a-h,o-z)            
      save
      parameter (tm32=2.d0**(-32),lenr=9689,ifdb=471)
      common/ransrb/ irs(lenr),next(lenr),ipoint,ipoinf
      irn=ieor(irs(ipoint),irs(ipoinf))   
      irs(ipoint)=irn                    
      rn=irn*tm32+0.7d0
      ipoint=next(ipoint)             
      ipoinf=next(ipoinf)            
      return                       
      end                         
